!> 4 阶 Runge-Kutta 法
module seakeeping_runge_kutta

    use seakeeping_kinds, only: wp => sk_real_kind
    implicit none

    private
    public :: runge_kutta

    abstract interface
        !> `fcn` 计算 ODE 的导数
        subroutine fcn(t, y, yp)
            import wp
            real(wp), intent(in) :: t       !! 自变量
            real(wp), intent(in) :: y(:)    !! 因变量
            real(wp), intent(out) :: yp(:)  !! 导数
        end subroutine fcn
    end interface

contains

    !> 4 阶 Runge-Kutta 法
    !> @note
    !> \[y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k4)\]
    !> 其中，\[k_1 = f(t_n, y_n)\]
    !> \[k_2 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1)\]
    !> \[k_3 = f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2)\]
    !> \[k_4 = f(t_n + h, y_n + hk_3)\]
    subroutine runge_kutta(func, t, y, dt, m, n)
        procedure(fcn) :: func                      !! 回调函数
        real(wp), intent(inout), dimension(*) :: y  !! 因变量
        real(wp), intent(inout) :: t                !! 自变量
        integer, intent(in) :: m                    !! 时间步数
        real(wp), intent(in) :: dt                  !! 时间步长
        integer, intent(in) :: n                    !! 因变量维数
        real(wp) :: h2
        real(wp), dimension(n) :: k1, k2, k3, k4
        integer :: i

        h2 = dt/2
        do i = 1, m
            call func(t, y(:n), k1)
            call func(t + h2, y(:n) + h2*k1, k2)
            call func(t + h2, y(:n) + h2*k2, k3)
            call func(t + dt, y(:n) + dt*k3, k4)
            y(:n) = y(:n) + dt*(k1 + 2*k2 + 2*k3 + k4)/6
            t = t + dt
        end do

    end subroutine runge_kutta

end module seakeeping_runge_kutta
